Superfluid density and phase relcixation in superconductors with strong disorder 
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As a prototype of a disordered superconductor we consider the attractive Hubbard model with 
on-site disorder. We solve the Bogoljubov-de-Gennes equations on two-dimensional finite clusters at 
zero temperature and evaluate the electromagnetic response to a vector potential. We find that the 
standard decoupling between transverse and longitudinal response does not apply in the presence 
of disorder. Moreover the superfluid density is strongly reduced by the relaxation of the phase 
of the order parameter already at mean-field level when disorder is large. We also find that the 
anharmonicity of the phase fluctuations is strongly enhanced by disorder. Beyond mean-field, this 
provides an enhancement of quantum fluctuations inducing a zero-temperature transition to a non- 
superconducting phase of disordered preformed pairs. Finally, the connection of our findings with 
the glassy physics for extreme dirty superconductors is discussed. 
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111 the last few years renewed interest emerged in the 
behavior of disordered superconductors at the verge of 
the metal-insulator transition. On the experimental side 
recent studies of disordered superconducting films at low- 
temperatures [1, 2] have revealed the existence of a strik- 
ing "pseudogap" behavior at large disorder, with the tun- 
neling conductance at zero bias starting to develop a sup- 
pression at temperatures much larger than T^. These 
findings suggest quite generically a separation between 
the energy scales associated to local pairing (gap and 
pseudogap) and superconducting phase coherence (super- 
fiuid density) with different dependencies on disorder. In 
particular, it is this second energy scale (to be more spe- 
cific, the superfluid stiffness) which is expected to control 
the stability of the superconducting phase in the proxim- 
ity to the superconductor-insulator (SC-I) transition. 

From the theoretical point of view the study of dis- 
ordered superconductors near the SC-I transition[3] has 
been based either on a bosonic approach[4], where the 
role of phase fluctuations emerges naturally, or on a more 
microscopic fermionic approach[5-8], which has put the 
focus on the emergence of a short-scale inhomogeneity 
induced by the strong disorder. In particular, in the lat- 
ter case the appearance of two characteristics features 
of the disordered SC state has been demonstrated: the 
spontaneous emergence of spatial structures in the local 
pairing gap and a general suppression of the phase coher- 
ence, which reflects in a decrease of the global superfluid 
stiffness[5-8]. These results appear already at mean-field 
level, as it has been shown by the Bogoljubov-de-Gennes 
(BdG) solution for a 2D SC system in the presence of 
on-site disorder [5], which serves also as comparison for 
more refined Monte Carlo results including thermal[7] or 
quantuiii[8] phase fluctuations. Recently it has been ar- 
gued that eventually a non self-averaging character typi- 
cal of glassy physics will appear at very strong disorder [9]. 



In the presence of disorder a non-trivial problem is 
posed by the correct computation of the superfluid stiff- 
ness Dg and by the understanding of the processes lead- 
ing to the supression of Dg. In the clean translational 
invariant case the superfluid behavior of the system re- 
flects the "rigidity" of the phase against an applied trans- 
verse vector potential. This leads to a purely dianiag- 
netic response of the current, which is the hallmark of 
the Meissner effect and the superfluid behavior of clean 
superconductors. At small disorder a paramagnetic re- 
sponse appears and the superfluid stiffness decreases with 
respect to the pure diamagnetic case. This paramagnetic 
current is the response of BdG quasiparticles to the ex- 
ternal vector potential, while keeping the phase of the 
order parameter unchanged. Within the BdG approach, 
one can then account for this effect by evaluating Dg as 
a disorder average over the BCS response function (the 
BCS current-current correlation function with no vertex 
corrections) , computed with the BdG solution for the dis- 
ordered system (cf. e.g. Ref. [5]). This procedure, which 
for small disorder is equivalent to the dirty-BCS limit, 
relies on the decoupling between longitudinal and trans- 
verse electromagnetic response: the phase of the order 
parameter does (does not) react to a longitudinal (trans- 
verse) field. However, this decoupling holds exactly only 
in clean systems. Our main results in this work are: i) 
For strong disorder we show that such a decoupling is 
strongly violated, leading to a dramatic decrease of the 
phase stiffness with respect to the dirty BCS case, due 
to the additional paramagnetic suppression coming from 
the phase relaxation to the applied transverse vector po- 
tential, a) We discuss the behavior of this paramag- 
netic phase response in connection to the formation of 
self-organized structures of the current in real space at 
scales much larger than the lattice spacing and the SC 
coherence length [cf. Fig. 1(b)], which share some analo- 



gies with the glassy features discussed in Refs. [9]. Hi) 
We show that the anharmonic phase fluctuations become 
strongly enhanced already at intermediate disorder with 
respect to the clean case, leading to a sizeable increase 
of quantum corrections to the superfluid stiffness due to 
phase fluctuations, with respect to a clean system. These 
results can account for the recently observed deviation of 
the superfluid stiffness from the dirty-BCS limit in NbN 
films[2], and offer new insight for the understanding of 
the SC-I transition. 

Our starting hamiltonian is the attractive Hubbard 
model with local disorder: 
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which we solve using the BdG equations [f 0] by allowing 
for a site dependent SC order parameter A^ — \U\{ci^Ci^) . 
The first sum is over nearest-neighbors pairs and we work 
on a i X L system using units such that lattice spacing 
and h ~ c ~ e ^ 1. The local potential Vi is randomly 
distributed between — Vq < Vi < Vq. We present re- 
sults for density n — 0.875 and at large SC coupling 
{U/t = —5), where the phase-relaxation effects and the 
connection to the glassy physics become more evident. 
Similar results are found in a wide parameter range. 

The superfluid stiffness Ds in the a — x,y direction 
(corresponding in the continuum limit to Dg — ns/m 
with Us being the superfluid density) is defined as the 
static limit of the transverse part, Kt, of the electromag- 
netic kernel, Kab, which describes the current response, 
ja ~ KabAh, to an applied electromagnetic field Ab (sum- 
mation over repeated indices is implicit) [11, 14]. 

In principle, the correct evaluation of the superfluid 
stiffness requires the knowledge of the electromagnetic 
kernel Kab within a so-called conserving approximation, 
i.e. an approximation that respects the gauge invariance 
of the theory. As it is well know, the BCS approximation 
for K^^^ is not conserving since, by neglecting the ver- 
tex corrections, it does not include the contribution of the 
phase relaxation to the current-current correlation func- 
tion [11, 14]. This is not a problem in the clean case as 
far as the transverse response is concerned, since phase 
fluctuations contribute only to the longitudinal part of 
-?i'a6[ll]- It is easy to check that this no longer holds in 
the presence of disorder. Indeed, upon considering also 
the non-local character of the BCS kernel for a disordered 
system, the effective action for phase fluctuations [16[ can 
in general be written at Gaussian level as: 
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In the absence of disorder the above expression simplifies 
considerably, since the BCS kernel K^^^ depends only 
on the difference r — r', so that by making the Fourier 
transform and approximating the BCS kernel with its 



(constant) long-wavelenght limit K^ 
action reduces to: 
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Integrating by parts the term V0 ■ A, one finds that 
couples only with the longitudinal part of the elec- 
tromagnetic field. Integrating out the phase one finds 
Ds = Kt = Kq'^^' i.e. the stiffness is given by the BCS 
kernel. Modeling disorder by a position dependent kernel 
in Eq. (3) one sees that the integration by parts leads to 
a coupling of the phase with both the longitudinal and 
the transverse part of the gauge field. The same is true 
if Eq. (2) is used which now depends separately on r and 
r'. For strong disorder one is not allowed to substitute 
the BCS kernel in Eq. (2) with its space-disordered aver- 
age which will restore translational invariance and erro- 
neously lead to the simplified Eq. (3). The consequences 
of this transverse- longitudinal decoupling breakdown are: 
(i) a change in the SC order-parameter phase even in the 
presence of a transverse field, contrary to the clean case; 
(ii) a failure of the BCS response function to compute 
the superfluid stiffness. 

In order to put these arguments on a quantitative basis 
we first show explicitly the behavior of the current in the 
presence of a constant vector potential A in the x direc- 
tion. Using the Peierls substitution, this corresponds to 
a change in the hopping term of Eq. (1) along x: 
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In a torus geometry (periodic boundary conditions) a 
constant A cannot be gauged away, and it corresponds 
to a flux <I> = AL through the torus. For a given dis- 
order configuration the microscopic current along x at 
finite A is then the sum of the paramagnetic and dia- 
magnetic contribution j^ = - %f^^ = (jf ) + (jf )• In 
the clean case jx can be directly derived from Eq. (3), 
and is proportional to the gauge-invariant phase gradi- 
ent W9 = V6 — 2A. Here 9 is the phase of the order 
parameter if we would eliminate A by a gauge transfor- 
mation. Since for an applied transverse field V0 — 0, 
then jj. reduces to the constant value —AD^. 

For the disordered system the result is radically differ- 
ent. If one neglects phase relaxation (BCS), each piece of 
the system responds according to its local stiffness, which 
is essentially determined by the local order parameter 
[Fig. l(a)[. Because this solution is not a saddle point, 
the current violates charge conservation and is clearly 
unphysical. Allowing for the phase to relax [Fig. l(b)[ 
both the current and the gauge invariant phase gradient 
[proportional to the line density in Fig. l(a)[ are strongly 
space dependent. This dependence is still correlated to 
the order parameter but now it is also conditioned by the 
existence of a percolative path. This means that the cur- 
rent response at one point can be strongly influenced by 
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FIG. 1: (Color online) Distribution of the local current (ar- 
rows) and lines of constant 9 superimposed to the map of the 
order parameter, neglecting phase relaxation (a) and allow- 
ing for it (b). We applied a transverse field A = 0.003, for 
a disorder amplitude Vo/t=2. Notice that in (b), while the 
percolative pattern of the current goes along the pattern of 
largest Ai values, the gauge-invariant phase gradient is larger 
along the locations of smaller Ai values. 



the local stiffness at very far points. Notice that isolated 
regions with a robust order parameter have practically 
zero current in panel (b). 

The gauge-invariant phase gradient V6' is large pre- 
dominantly in regions with minimum values of the order 
parameter A^. This can be understood by mapping the 
problem to a random resistor network[12, 13] where the 
bad SC regions map into poor conducting regions and 
9 maps into the electrostatic potential, so large "poten- 
tial" drops concentrate in the "poor" conducting regions. 
Moreover, in contrast to the clean isotropic case, where 
jx has the same direction than the gauge-invariant \/9 
(see Eq. (3)), here one expects a tensorial relation be- 
tween the two, as given in Eq. (2) and confirmed by the 
results of Fig. 1. It is worth stressing that while the local 
SC order parameter can be quite small, the local density 
of states always shows a sizable gap, as found in various 
previous analyses[5, 8], supporting the view that coher- 
ent superconductivity takes place in Fig. 1 on a system 
of localized preformed pairs. 

A direct procedure to compute the superfluid stiffness 
beyond BCS that accounts for the SC-phase relaxation, 
like those shown in Fig. 1, is given by the change in the 
ground-state energy E in the presence of the constant 
vector potential of Eq. (4), [14, 15] 
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The resulting values of Ds are shown in Fig. 2, along with 
their BCS counterparts. While the effects seem not dra- 
matic at the scale of the bare stiffness they are gigantic 
when the stiffness is small due to strong disorder, as can 
be seen in the right inset. Such a dramatic difference is 
clearly due to the global reorganization of the current, 
i.e. only good SC regions that are along the percolative 
path contributes to the stiffness. In the BCS case instead 
the global stiffness results to be a simple spatial average 
of the local stiffness. 



As the stiffness decreses due to increasing disorder also 
quantum phase fluctuations beyond Gaussian level be- 
come relevant [5, 16]. In general, in the presence of a 
constant gauge fleld A the expansion of the energy per 
unit surface, up to fourth order can be written as: 
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where 59 is the phase gradient in a given {x or y) direc- 
tion and we have restored explicitly the lattice spacing a. 
From Eq. (6) it is evident that while Ds provides the stiff- 
ness the coupling Q measures the anharmonicity of the 
phase fluctuations. The second part of Eq. (6) translates 
the expansion into an analogus expression for the phase 
gradient by means of the mininal-coupling substitution 
A -^ 69/2. 

In Eq. (6) we also introduced the couplings J = 1)^/4 
and J' = (3a^/(4^Q) in order to make a closer analogy 
with the XY model, H = JJ2{tj)i'^ - cos(6'j - 9j)) that 
is usually considered as the prototype model for phase 
fluctuations in a superconductor. The latter is deflned in 
a coarse-grained scale ^o given by the coherence length. 
Notice that contrary to J, the characteristic energy for 
anharmonic fluctuations J' in Eq. (6) depends on the 
coarse-grained scale ^o- For the XY model the gradient 
expansion is given by Eq. (6) with J' — J. 

Quantum corrections to D^ can be computed with a 
perturbative approximation from Eq. (6): 



D. = D^. 
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where D'^ is the (bare) value obtained by the Gaussian 
expansion in Eq. (6) and {{69?) is computed including 
the dynamics of the phase (cf. e.g. ]16, 17]). 

The coherence length has been determined from the de- 
cay of the gap amplitude correlations which yields ^o ~ 
la, so that for the present parameters Q/Dg ~ J' /J. The 
upper left inset to Fig. 2 shows the ratio Q/Dg as obtained 
from the fit of E{A) to Eq. (6). For increasing disorder 
Q/Dg « J' /J becomes strongly enhanced leading rapidly 
to a larger phase-fluctuation correction within the Hub- 
bard as compared to the XY model where J'/J = 1. 
Above a critical value of the disorder strength Vo/f w 3 
a SC-I transition is obtained when Dg = 0. Within an 
RPA analysis of the full model on small clusters we have 
also checked [18] that {{69?) is increasing with disorder, 
however, the dependence on Vq being much weaker than 
that for Q/Dg. 

The above discussion shows that the XY model fails 
in general to describe quantitatively corrections of the 
superfluid stiffness due to numerical differences in the 
quartic coefficient. Notwithstanding a disordered version 
of it can help to get some insight on the anomalously 
large increase of the J' / J ratio as a function of disor- 
der. Let us consider for simplicity a one-dimensional 
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FIG. 2: (Color online) Main panel: Comparison between the 
superfluid stiffness computed using Eq. (5) and its BCS coun- 
terpart as a function of disorder. Upper riglit inset: Same 
data plotted on a log-log scale to demonstrate increasing de- 
viations for large Vo/t. Upper left inset: Ratio between an- 
harmonic (Q) and quadratic coefficient (Ds) of the expansion 
Eq. (6) (with Q/Ds = 1.01 at V = 0). Note that for the 
given parameter set ^o ~ la so that J' /J ~ Q/Ds- 



XY-model in the presence of a constant gauge field, 
H = J2iJi,i+i^i^ - cos(6'j - 0i+i - 2A)). In order to 
describe a disordered model we consider a random dis- 
tribution of couplings Jis+i- To derive the coefficients 
J, J' of the expansion (6) we can use the fact that in 
ID the current in each link, j(v4.), must be conserved, so 
that —Ji^i+i sin(6',i — 6'i+i — 2A) = i{A)/2. By inverting 
this relation, summing over the index i and expanding 
j{A)/2 « J A ~ ^A'^ we then obtain that the coefficients 
are given by the spatial averages, J^^ 

J4 



l/Ji^i+i (see 

also Ref. [12]), and J' — J* '^/Jfi+i- By computing 
for example the ratio J' / J for a Gaussian distribution of 
local coupling values one can see that it increases very 
rapidly as the distribution width increases, leading to a 
non- vanishing occurrence probability for very low values 
of Ji,i+i, which have a much stronger effect on J' than 
on J . As a byproduct this computation also shows the 
importance of phase relaxation, without which J would 
be defined as J = Ji^i+i- Obviously, in 2D the current in 
each link is not constant. Thus in principle the current 
would try to avoid links with very low J.y couplings to 
avoid coexistence of large phase change and large current, 
which from Fig. 1 seems to be indeed the case. However, 
these patterns have an almost one dimensional character 
so that the above argument could be still qualitatively 
correct even in 2D. Nevertheless the variation of local Jij 
values along the active paths (i.e. paths with current) 
can be smaller than along paths with vanishing current. 
In conclusion, we have analyzed from a real space BdG 
approach for the attractive Hubbard model, the reaction 
of a strongly disordered SC system to an applied trans- 
verse vector potential. The first main result concerns the 
standard decoupling between transverse and longitudinal 



response which does not apply in the presence of disorder. 
We have evaluated the corrections of the superfluid stiff- 
ness to the dirty BCS limit (cf. Fig. 2) and obtained a 
strong reduction of the superfluid stiffness due to the re- 
laxation of the order parameter phase for large disorder. 
We also find that the anharmonicity of the phase fluctu- 
ations is increased by disorder. Beyond mean-field, this 
enhances the effect of quantum fluctuations in producing 
a zero-temperature transition to a non-SC phase. Both 
effects can explain the deviations of the zero-temperature 
superfluid stiffness from the BCS dirty limit reported re- 
cently for strongly disordered NbN films [2]. 

Secondly, our calculations have revealed that the stiff- 
ness for strong disorder is dominated by quasi one- 
dimensional percolative paths along sites with large gap 
parameters. In some sense this mimics the finding of 
Ref. [9] where the phase diagram of a model, similar 
to the present one in the large U /t limit, has been ana- 
lyzed within the so-called cavity mean-field approxima- 
tion. Interestingly it was found that there is a regime 
of broken-replica symmetry where the partition function 
is determined by a small number of paths. Taking that 
a similar glassy physics is at work in our model, where 
the response is governed by a few number of percolative 
paths, then this should also appear in the distribution 
functions of order parameters and energy level spacings. 
An analysis of this issue is in progress. 
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